Shallow water dynamics

This notebook describes the shallow water model and phenomena that can be investigated with it.

Shallow water equations

The following relies heavily on material in Atmospheric and Oceanic Fluid Dynamics by Geoffrey K. Vallis (AOFD).

The shallow water equations describe a fluid layer of constant density in hydrostatic balance, bounded by a rigid surface at the bottom and a free surface at the top.

The equations for conservation of mass $(1)$ and momentum $(2)$ are (note that $\frac{d}{dt} = \frac{\partial}{\partial t} + \pmb{v} \cdot \nabla$):

\begin{equation} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \pmb{v}) = 0 \label{eq: 1}\tag{1} \end{equation}\begin{equation} \frac{d\pmb{v}}{dt} + 2\pmb{\Omega} \times \pmb{v} = -\frac{1}{\rho}\nabla p - \nabla \Phi \label{eq: 2}\tag{2} \end{equation}

Momentum

The z-component of $(2)$, using $\Phi = gz$ and assuming hydrostatic equilibrium is:

\begin{equation} \frac{\partial p}{\partial z} = -\rho g \label{eq: 3}\tag{3} \end{equation}

As $\rho$ is constant, we can integrate this. If we use the boundary condition that $p = 0$ at the surface of the fluid, $h_s$ i.e. when $z = h_s(x, y, t)$, we get:

\begin{equation} p(x, y, z, t) = \rho g (h_s(x, y, t) - z) \label{eq: 4}\tag{4} \end{equation}

Let us take $\pmb{v} = (u, v, w)$ and $\pmb{\Omega} = \Omega (0, cos(\theta), sin(\theta))$, where $\Omega$ is the rotation rate of the earth and $\theta$ is latitude. Then the x and y components of $(2)$, using $(4)$ and neglecting the coriolis term containing vertical velocity, $w$, become:

\begin{equation} \frac{du}{dt} - fv = -g\frac{\partial h_s}{\partial x} \\ \frac{dv}{dt} + fu = -g\frac{\partial h_s}{\partial y} \label{eq: 5}\tag{5} \end{equation}

where $f = 2\Omega sin(\theta)$. In the model, we consider a beta plane, so $f = f_0 + \beta y$ where $f_0 = 2\Omega sin(\theta_0)$ and $\beta = (df/dy)_{\theta_0} = 2\Omega cos(\theta_0) / R_{planet}$ using the approximation $y = R_{planet}\theta$.

Mass conservation

As $\rho$ is constant, $(1)$ becomes:

\begin{equation} \nabla \cdot \pmb{v} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0 \label{eq: 6}\tag{6} \end{equation}

From $(5)$, $u$ and $v$ are independent of $z$ as nothing in $(5)$ depends on $z$. We can thus integrate $(6)$ over $z$ from the bottom rigid surface, $z = h_b(x,y)$ to the top surface, $z = h_s(x,y,t)$, giving:

\begin{equation} w(h_s) - w(h_b) + (h_s - h_b)\left[\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right] = 0 \label{eq: 7}\tag{7} \end{equation}

The vertical velocity is just the material derivative of the position of a fluid element so $w(h_s) = \frac{dh_s}{dt}$ and $w(h_b) = \frac{dh_b}{dt}$. So if we define the fluid depth, $h = h_s - h_b$, we have:

\begin{equation} \frac{dh}{dt} + h\left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right) = 0 \label{eq: 8}\tag{8} \end{equation}

Conservative form

For the model, we use equations $(5)$ and $(8)$ in the conservative form shown below (can be found through algebra on $(5)$ and $(8)$):

\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial (uh)}{\partial x} + \frac{\partial (vh)}{\partial y} = 0 \\ \frac{\partial (uh)}{\partial t} + \frac{\partial (u^2h + gh^2/2)}{\partial x} + \frac{\partial (uvh)}{\partial y} = h\left(fv - g\frac{\partial h_b}{\partial x}\right) \\ \frac{\partial (vh)}{\partial t} + \frac{\partial (uvh)}{\partial x} + \frac{\partial (v^2h + gh^2/2)}{\partial y} = h\left(-fu - g\frac{\partial h_b}{\partial y}\right) \\ \label{eq: 9}\tag{9} \end{equation}

This can be written in vector form as:

\begin{equation} \frac{\partial \pmb{U}}{\partial t} + \frac{\partial F(\pmb{U})}{\partial x} + \frac{\partial G(\pmb{U})}{\partial y} = Q(\pmb{U}) \label{eq: 10}\tag{10} \end{equation}

Where:

\begin{equation} \pmb{U} = \begin{bmatrix} h \\ uh \\ vh \end{bmatrix} , \; F(\pmb{U}) = \begin{bmatrix} \pmb{U}_1 \\ \pmb{U}_1^2/\pmb{U}_0 + g\pmb{U}_0^2/2 \\ \pmb{U}_1 \pmb{U}_2 / \pmb{U}_0 \end{bmatrix} , \; G(\pmb{U}) = \begin{bmatrix} \pmb{U}_2 \\ \pmb{U}_1 \pmb{U}_2 / \pmb{U}_0 \\ \pmb{U}_2^2/\pmb{U}_0 + g\pmb{U}_0^2/2 \end{bmatrix} , \; Q(\pmb{U}) = \begin{bmatrix} 0 \\ \pmb{U}_0 \left(f\pmb{U}_2 / \pmb{U}_0 - g\frac{\partial h_b}{\partial x}\right) \\ \pmb{U}_0 \left(-f\pmb{U}_1 / \pmb{U}_0 - g\frac{\partial h_b}{\partial y}\right) \end{bmatrix} \label{eq: 11}\tag{11} \end{equation}

The model solves the equations by putting the physical values in this vector form and then poviding the relavent flux and source functions to a chosen numerical solver.

The numerical solvers are in the file numerical_methods.py and are based on this video. The default method is richtmyer.

Example simulations

Various phenomena can be visualised by using the function run_simulation defined in the cell below. This function has been initialised with default values which tend to give quite good results. The only non-default argument is initial_info. This is a dictionary which specifies the initial values of $h_s$, $u$ and $v$ for the simulation. There are four options for the type of simulation, indicated by initial_info['type']. These are 'uniform_zonal', 'sinusoidal_zonal', 'jet_zonal', 'height_gaussian' or 'height_step'. Each of these requires different information as described here.

The phenomena displayed here are inspired by the computer practical by Paul Connolly and Robin Hogan. Each requires different initial_info and some are better displayed with parameters which differ from the default.

Geostrophic adjustment

This example illustrates the phenomena of geostrophic adjustment described in chapter 3.9 of AOFD.

Geostrophic adjustment is when a fluid in an initially unbalanced state, naturally evolves to a state of geostrophic balance. To show this, we consider a fluid initially at rest but with a discontinuity in the height field at x=0.

No rotation

The cell below shows the case with no rotation. You can see that the disturbance is radiated away to infinity at speed $c = \sqrt{gh_{mean}}$. So with $x < ct$, $h_s = h_{mean}$ and $u = -(h_{s, max}-h_{mean})g/c$. With $x > ct$, $h_s = h_{s, initial}$ and $u = 0$.

With rotation

The next cell shows the case with rotation. In this case, the disturbance is not radiated away to infinity. Potential vorticity adjustment constrains the adjustment to within a deformation radius, $L_d = c/f_0$. Within the deformation radius (in this example a deformation radius is the distance between a pair of adjacent velocity arrows), $h_s$ goes through zero and the veritcal velocity $v$ is positive. Far away from the initial disturbance, the fields retain their initial values at all times.

Now the bottom plot shows non-zero relative vorticity, $\zeta$ too. In particular, as potential vorticity, $q = (\zeta + f)/h$ is conserved, where h is lower than it started, just to the right of the initial disturbance, $\zeta$ must drop to keep $q$ constant and vice-versa to the left.

Increased deformation radius

The below cell shows a simulation with a deformation radius four times as large. At times not affected by the waves such as at 1 day, you can see that $v$ is positive over a greater region, the transition of $h$ from low to high is smoother and the vorticity is non zero over a greater area.

Gravity waves

If we turn rotation off and introduce a gaussian blob onto the surface height field, we can illustrate gravity waves which propogate away from the initial disturbance at speed $c = \sqrt{gh_{mean}}$.

With orography

If we add some structure to the floor we can see what happens for example in a Tsunami, where a wave enters a region of shallower water. Here we introduce a gaussian mountain on the floor, indicated by the green contours.

You can see that the height high and low extremes are larger around the shallow region in the centre. You can also see the waves slow down as they enter the shallower water causing a disfiguration of the circular annulus.

Barotropic instability

Following chapter 9.3 of AOFD, barotropic instability is an instability that arises in a fluid with a shear in the flow i.e. if the x velocity, $u = u(y)$, is a function of only $y$.

The condition for instability can be obtained by subbing in $u = u(y) + u'(x, y, t)$, $v = v'(x, y, t)$ and $h = h(x, y) + h'(x, y, t)$ into $(5)$ and $(8)$ where the dashed quantities are small perturbations. $(8)$ then gives $\frac{\partial u'}{\partial x} + \frac{\partial v'}{\partial y} = 0$ so we can define $u'$ and $v'$ from a stream function: $u' = -\frac{\partial \psi '}{\partial y}$ and $v' = \frac{\partial \psi '}{\partial x}$. Subbing this into $(5)$, we can get a single equation for $\psi '$. Then seeking a solution of the form $\psi ' = \phi (y) e^{ik(x-ct)}$, instability arises when the imaginary part of $c$ is positive.

It can be found that a necessary condition for this instability is that gradient of absolute vorticity, $\zeta_a = f + \zeta_{mean}$, changes sign at least once in the domain ($\zeta_{mean}$ is the relative vorticity of the mean flow, so in this case with $u = u(y)$, $v=0$ and $f = f_0 + \beta y$: $\zeta_a = \beta - \frac{\partial u}{\partial y}$).

Jet

To simulate this, we start the simulation by adding random noise to a height field whch is such that in geostrophic equilibrium, there is a zonal jet at y=0. It is clear from the vorticity plot at $t=0$ that the gradient changes sign near the jet so should lead to instability.

Note that without adding the random noise, nothing happens because the initial condition is a steady solution to the equations of motion.

A good explanation of barotropic instability is given in chapter 5.4 of Atmospheric Circulation Dynamics and General Circulation Models in terms of an interaction between Rossby waves.

Jupiter

In this example, we start with zonal structure such that the zonal wind oscillates between westerly and easterly as on Jupiter. The initial height field is then computed such that it is in geostrophic equilibrium.

Adding noise to this and then running the simulation leads to the development of large structures, which look somewhat like the giant red spot on Jupiter. This arises because (as explained in chapter 11.3 of AOFD) in 2D turbulence, small scale sructures merge into larger ones. There is an upwards cascade of energy.

Mountain waves

The next example illustrates the effect of a mountain on a uniform westerly flow with the height field set initially to be in geostrophic balance.

Gravity

Initially, we set $\beta = 0$ so just consider the case where the restoring force is provided by gravity. We can see in the cell below a meandering that occurs in the lee of the mountain. This is because the fluid depth, $h$ is lower over the mountain so potential vorticity, $q = (f + \zeta)/h$ increases. To offset this, the relative vorticity, $\zeta$ must decrease. This causes the anticylonic motion just beyond the mountain. Then as the fluid moves beyond the mountain, $h$ starts to increase and thus $q$ starts to decrease again. To compensate for this, $\zeta$ must increase. This then pushes the fluid into a region where $h$ is lower and the process repeats.

Rossby

If we set $\beta \neq 0$, then we also include the restoring force due to the coriolis effect. Now when the fluid moves to negative y just beyond the mountain, $f$ decreases as well as $h$ increasing. These effects both act to decrease $q$ so $\zeta$ must increase more rapidly to compnensate. The effects of $f$ and $h$ always act in phase so $\zeta$ must now change more sharply, resulting in more meanders.

Equator

At the equator, we substitute $f = \beta y$ into $(5)$ and the process of deriving the waves is explained in chapter 4.7.3 of Atmospheric Circulation Dynamics and General Circulation Models. It is shown that equatorial waves decay exponentially away from the equator with the fall off being set by the deformation radius $L_d = \sqrt{c/(2\beta )}$ where $c = \sqrt{gh_{mean}}$.

The first plot below shows that for a Jupiter like experiment performed at the equator, the instabilities that form outside the equator do not enter the equator.

Kelvin and Rossby Waves

The next simulation shows that a disturbance centered at the origin produces an eastward propogating Kelvin wave and a westward propogating Rossby wave. The Rossby wave moves slower and the height field peaks off the equator. When the Kelvin wave hits the wall, you can see that it is reflected back as Rossby waves and paritially propogates poelward as Kelvin coastal waves.

A similar simulation is shown in AOFD figure 22.18.